{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!--BOOK_INFORMATION-->\n",
    "<img align=\"left\" style=\"padding-right:10px;\" src=\"figures/PHydro-cover-small.png\">\n",
    "*This is the Jupyter notebook version of the [Python in Hydrology](http://www.greenteapress.com/pythonhydro/pythonhydro.html) by Sat Kumar Tomer.*\n",
    "*Source code is available at [code.google.com](https://code.google.com/archive/p/python-in-hydrology/source).*\n",
    "\n",
    "*The book is available under the [GNU Free Documentation License](http://www.gnu.org/copyleft/fdl.html). If you have comments, corrections or suggestions, please send email to satkumartomer@gmail.com.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!--NAVIGATION-->\n",
    "< [Evaporation](04.05-Evaporation.ipynb) | [Contents](Index.ipynb) | [Surface Water](04.07-Surface Water.ipynb) >"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.6 渗透量\n",
    "\n",
    "Green-Ampt法给出的累计渗透量写为:\n",
    "\n",
    "<center>$F(t)-\\psi\\Delta\\theta ln(1+\\frac{F(t)}{\\psi\\Delta\\theta})=Kt\\quad\\quad\\quad\\quad(4.12)$</center>\n",
    "\n",
    "其中，$F(t)$是t时刻后的累计渗透量，$\\psi$是吸头，$\\Delta\\theta$给出为，\n",
    "\n",
    "<center>$\\delta\\theta = (1-S_e)\\psi_e\\quad\\quad\\quad\\quad\\quad\\quad\\quad\\quad\\quad(4.13)$</center>\n",
    "\n",
    "其中，$S_e$是饱和度，$\\theta_e$是有效的孔隙度，$K$是渗透系数。用迭代过程求解方程，方程4.12被重写为，\n",
    "\n",
    "<center>$F(t) = \\psi\\delta\\theta ln(1+\\frac{F(t)}{\\psi\\delta\\theta})\\quad\\quad\\quad\\quad\\quad\\quad\\quad(4.14)$</center>\n",
    "\n",
    "我们使用`while`函数来迭代，直到我们达到了求的精度。使用`append`方法来存储$F$的迭代值。`append`通过一个个的项追加数组，并将输入变量放入其中。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import division\n",
    "import numpy as np\n",
    "# define the variables\n",
    "theta_e = 0.486\n",
    "psi = 16.7\n",
    "K = 0.65\n",
    "S_e = 0.3\n",
    "t = 1\n",
    "\n",
    "#calculate dtheta\n",
    "dtheta = (1-S_e)*theta_e\n",
    "\n",
    "# initial guess of F\n",
    "F_old = K*t\n",
    "epsilon = 1\n",
    "F = []\n",
    "while epsilon > 1e-4:\n",
    "    F_new = psi*dtheta * np.log(1+F_old/(psi*dtheta)) + K*t\n",
    "    epsilon = F_new - F_old\n",
    "    F_old = F_new\n",
    "    F.append(F_new)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在，我们绘制$F$的迭代值的图，看看$F$如何用迭代来更新。我们也在`plot`函数中使用`-ok`。`-o`代表示填充点的连线，`k`表示图的颜色为黑色。我们还指定了`xlabel`和`ylabel`的字体大小。我们已经为`ylabel`使用'25'，为`xlabel`使用’20‘来表示不同的字体大小可以用于不同的文本。当然，没有必要为`ylabel`和`xlabel`定义不同的字体大小。同样的参数`fontsize`也可以用来定义图例的字体大小。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEYCAYAAACUdWs9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xl8VOXd///XJ2GN4AqiN5LEW6ytG4gBqtgKVSta64K9FX4pd9XaVLu5td5YqqVavG1tXWvFuGF/Rr1rjRYVrVoTt4osiohiVZRNUFBWAYWEz/ePc0aHycxkm8yZmbyfj8c8ZuY613XOZ85Mzifnus5i7o6IiEgmFEUdgIiIFA4lFRERyRglFRERyRglFRERyRglFRERyRglFRERyRglFRERyZhIk4qZ9TCzmWb2qpm9bma/SVLnQjN7w8zmmdk/zawsblqjmc0NH9OyG72IiCSyKE9+NDMDdnD3T8ysK/A8cJ67z4irMwp4yd03mdm5wEh3Pz2c9om794okeBERaaJLlAv3IKN9Er7tGj48oU5d3NsZwHfbs8w+ffp4eXl5e2YhItKpzJkz5yN379uSupEmFQAzKwbmAAOBm9z9pTTVvw88Fve+h5nNBhqAq9z9oeaWV15ezuzZs9sTsohIp2Jmi1taN/Kk4u6NwGAz2xl40MwOdPf5ifXM7LtABXBkXHGpuy83s/8Enjaz19x9YZK2VUAVQGlpaYd8DhERyaGjv9x9LVAPjE6cZmZHAxOBE939s7g2y8Pnd8O2h6SYd7W7V7h7Rd++LdqDExGRNoj66K++4R4KZtYTOBp4M6HOIcAtBAllZVz5LmbWPXzdBxgBvJGt2EVEpKmou7/2BO4Kx1WKgL+6+yNmdjkw292nAVcDvYD7g4PFWOLuJwJfAW4xs21h26vcXUlFRCRCUR/9NY8kXVbuflnc66NTtP0XcFDHRSciIq2VM2MqIiIxNTU1lJeXU1RURHl5OTU1NQXTLh9ibBd371SPQw891EU6s7vvvtvLysrczLysrMzvvvvuDmnTnmWVlJQ4wTlrDnhJSUmzbfOhXT7EmAzBcESLtrGRnlEfhYqKCtd5KpJLampqmDhxIkuWLKG0tJTJkydTWVnZIe1qamqoqqpi06ZNn5eVlJRQXV2dsm26NuPGjaOhoeHzR2Nj4+ev77//fiZMmMDmzZs/b9ejRw8uvfRSRo8ezbZt22hsbGzyPHbsWFauXNkkjr59+3LLLbfg7mzbtu3z59jrn/3sZ3z88cdN2u26665cddVVKTeCl156KWvWrGnSbpddduHSSy8lto2Mf3Z3rrzyStauXduk3c4778zFF1/cpA3A1Vdfzbp165q02WmnnbjoootI3B7H3l977bUp251//vlN2sXa3nDDDUnblZWVsWjRoiblqZjZHHevaFFdJRWR6LRlIw9w9913U1VVtd0Gu3v37lx44YUcdthhbN68+fPHpk2bPn99ww03sGHDhibz69GjB4cffjhbtmxp8njvvfdobGzM7AeXSJkZ27Zta019JZVUlFSko7Rmz8Hd2bhxI/vttx/Lly9vMr13796cfvrprF+/nnXr1m33HHu0VlFRUdoNyRFHHEG3bt2aPO65556UbX7zm9/QpUsXiouL6dKly3aPc845J2kbM+Ohhx6iqKiIoqIiiouLt3seO3YsH374YZN2e+65J9OnT6eoqAgz+7x97PWoUaOSrsv+/fszY8YMzCzp49BDD2XZsmVN2u2111689tprn9eLxR573n///Vm6dGmTdgMGDOCtt95K2mbgwIEsWbKkSZvS0lLefffd7erHr6/y8vKk7Zrb4ygvL2fx4qYnw3fknkrkYxzZfmhMRZqTqXGAbt26eWVlpU+YMMHPOussP+GEE3zo0KFeVlbmPXv23K5usseee+7p++23nw8bNsyPPvpoHzNmjJ955pl+3nnnpWxjZj5z5kyfP3++L1y40JcvX+5r1qzxzz77zLdt2+ZlZWVJ25WVlaX8bG1p0552+TLuoDGV5I/IN/LZfiipSDrN/RE2NDT44sWL/emnn/bbbrvNL7nkEj/ttNO8W7duKTf0Xbp08f/4j//wwYMH+7HHHuvjx4/3iy66yH//+9/7brvtlvMb7Cg2aNk8MCDb7fIhxkRKKkoq0kapNtY9evTwL33pS02SR5cuXXyfffZJu+ewbdu2lMvLlw12lBs0iZ6SipKKeMs3aB999JE/9thj/pvf/CZtd9Spp57qF198sd9yyy3+5JNP+rvvvutbt25197bvObQmzky1E2mt1iQVDdRLQUp1VNWNN97Ivvvuy6xZs5g5cyazZs3aboC0uLiYhoaGJvNrbmCzrUdxieQDHf2VhpJK55DqqJd4paWlDB06lGHDhjF06FAOPfRQHn744TYnh7aebyKS65RU0lBSKXwLFy5k4MCBKac//PDDDB06lH79+iWdruQgsj0llTSUVAqPu/P6669TW1vLAw88wLx581LWbe3x+SLSuqSiC0pKXkh2UbxZs2ZxySWX8OUvf5mDDjqISZMmseOOO3LNNddw3XXXUVJSst08SkpKmDx5ckSfQKRziPp+KiLNShwEX7x4MePHj8fdKS4u5hvf+AYXXHABJ510Envuuefn7fr06aNuLJEsU/eX5LzS0tKkl8PYbbfdeOutt9h1110jiEqk82hN95f2VCRnbdiwgT//+c9JEwrA6tWrlVBEckzU96jvYWYzzexVM3vdzH6TpE53M/s/M3vHzF4ys/K4aZeE5f82s2OzGbt0nPXr13PllVdSXl7OhAkT6NGjR9J6paWlWY5MRJoT9UD9Z8A33H0QMBgYbWZfTajzfWCNuw8ErgV+B2Bm+wNjgQOA0cCfw3vdS55au3Ytl19+OWVlZUycOJHDDjuMGTNmcNttt2nQXSRPRJpUwisAfBK+7Ro+Egd5TgLuCl//DTjKgmtDnwTc5+6fuft7wDvAsCyELRm2evVqLrvsMsrKyvj1r3/NkUceyezZs3nkkUcYPnw4lZWVVFdXU1ZWhplRVlamM9VFclTkYyrh3sUcYCBwk7u/lFClP7AUwN0bzGwdsFtYPiOu3rKwTHJY/ImF/fv3Z8iQIdTV1bFhwwbGjBnDpZdeyuDBg5u0q6ysVBIRyQORJxV3bwQGm9nOwINmdqC7z4+rYsmapSlvwsyqgCpQP3yUEg8NXrZsGcuWLWP48OHceuutHHTQQRFHKCLtFfWYyufcfS1QTzA+Em8ZMADAzLoAOwGr48tDewFNb/sWzLva3SvcvaJv374ZjlxaauLEidtdUyvmgw8+UEIRKRBRH/3VN9xDwcx6AkcDbyZUmwZ8L3z9HeDp8FLM04Cx4dFhewP7AjOzE7m01meffZbyAo/JbpMqIvkp6u6vPYG7wnGVIuCv7v6ImV1OcP3+acDtwP9vZu8Q7KGMBXD3183sr8AbQAPw47ArTXLMwoULOe2001JOV5ekSOGINKm4+zzgkCTll8W9/hT4rxTtJwM6rjSH3X///Zx99tkUFRVxwQUXcMsttzS5rLwODRYpHDkzpiKF5dNPP+VHP/oRp512Gl/5yleYO3cu11xzjQ4NFilwuvaXZNzbb7/Naaedxty5c7nooou48sor6datW9RhiUgb6dpfEpn77ruPH/zgB3Tr1o2HH36YE044IeqQRCSL1P0lGbF582Z++MMfMm7cOA4++GDmzp2rhCLSCSmpSJvE3zSrf//+7LffflRXV3PxxRdTX1/PgAEDmp+JiBQcdX9JqyWeGb98eXDO6c9//nN+97vfRRmaiERMeyrSaqnOjL///vsjiEZEcomSirRaqjPgdWa8iCipSKvtsssuSct1ZryIKKlIqzz00EOsXr2a4uLt74emM+NFBJRUpBVeeOEFxo0bx7Bhw3RmvIgkpaO/pEXeeOMNvv3tb1NaWsqjjz5Knz59OOuss6IOS0RyjPZUpFnvv/8+o0ePpnv37jz++OP06dMn6pBEJEdpT0XSWrt2LaNHj2bt2rU888wz7L333lGHJCI5TElFUvr00085+eST+fe//81jjz3GIYc0uUuBiMh2lFQkqcbGRsaPH88zzzzDPffcw1FHHRV1SCKSBzSmIk24OxdccAF/+9vf+OMf/8i4ceOiDklE8kSkeypmNgD4C7AHsA2odvfrE+r8Aogdq9oF+ArQ191Xm9kiYAPQCDS09Hr/kt7vf/97brzxRi688EIuvPDCqMMRkTwSdfdXA3CRu79sZr2BOWb2pLu/Eavg7lcDVwOY2beBC9x9ddw8Rrn7R1mNuoD95S9/YcKECYwbN46rr7466nBEJM9E2v3l7ivc/eXw9QZgAdA/TZNxwL3ZiK2ziL+Efb9+/TjjjDM46qijmDp1KkVF6h0VkdbJma2GmZUDhwAvpZheAowGHogrduAJM5tjZlUdHWOhiV3CfvHixbg7K1euBGDs2LG6/a+ItElOJBUz60WQLM539/Upqn0beCGh62uEuw8BjgN+bGZfTzH/KjObbWazV61aldHY81myS9i7O7/97W8jikhE8l3kScXMuhIklBp3r01TdSwJXV/uvjx8Xgk8CAxL1tDdq929wt0r+vbtm5nAC4AuYS8imRZpUjEzA24HFrj7NWnq7QQcCfw9rmyHcHAfM9sB+CYwv2MjLiypLlWvS9iLSFtFvacyAhgPfMPM5oaP483sHDM7J67eKcAT7r4xrqwf8LyZvQrMBB5198ezF3r+u/jii5uU6RL2ItIekR5S7O7PA9aCelOBqQll7wKDOiSwTuLFF1+kuLiYfv36sWLFCkpLS5k8ebIuYS8ibRb1eSoSkbq6Ou6++25+9atfccUVV0QdjogUCHP3qGPIqoqKCp89e3bUYURqy5YtDBo0iC1btjB//nx69uwZdUgiksPMbE5Lr1iiPZVO6A9/+ANvvvkm06dPV0IRkYyKeqBesuy9997jiiuuYMyYMRx33HFRhyMiBUZJpRNxd372s59RXFzMddddF3U4IlKA1P3Vifz973/nkUce4Q9/+AMDBgyIOhwRKUAaqO8kPvnkE/bff3923nln5syZQ9euXaMOSUTyhAbqpYnLL7+cpUuXcu+99yqhiEiH0ZhKJzB//nyuvfZavv/97zNixIiowxGRAqakUuC2bdvGueeey0477cTvfve7qMMRkQKn7q8Cd9ddd/H8889z++23s9tuu0UdjogUOO2pFLCPP/6YX/ziFxx++OGcccYZUYcjIp2AkkoBmzBhAmvXruXmm2/WrYFFJCu0pSlQ//rXv7jttts4//zzOfjgg6MOR0Q6CSWVAtTQ0MC5557LXnvtxaRJk6IOR0Q6EQ3UF6AbbriBefPm8cADD9CrV6+owxGRTkR7KgWipqaG8vJyioqK+PnPf86gQYM45ZRTog5LRDqZqO9RP8DM6sxsgZm9bmbnJakz0szWxd1u+LK4aaPN7N9m9o6ZTchu9LmjpqaGqqoqFi9ejLvj7rz11lvcc889UYcmIp1MpNf+MrM9gT3d/WUz6w3MAU529zfi6owEfu7uJyS0LQbeAo4BlgGzgHHxbZMpxGt/lZeXs3jx4iblZWVlLFq0KPsBiUhBac21vyLdU3H3Fe7+cvh6A7AA6N/C5sOAd9z9XXffAtwHnNQxkea2JUuWtKpcRKSj5MyYipmVA4cALyWZfJiZvWpmj5nZAWFZf2BpXJ1ltDwhFZTS0tJWlYuIdJScSCpm1gt4ADjf3dcnTH4ZKHP3QcCNwEOxZklmlbQvz8yqzGy2mc1etWpVpsLOGZMnT6a4uHi7spKSEiZPnhxRRCLSWUWeVMysK0FCqXH32sTp7r7e3T8JX08HuppZH4I9k/g7Te0FLE+2DHevdvcKd6/o27dvxj9D1I499ljMjF69emFmlJWVUV1dTWVlZdShiUgnE+l5KmZmwO3AAne/JkWdPYAP3d3NbBhBIvwYWAvsa2Z7A+8DY4H/LzuR55Zbb72VhoYGZsyYwQEHHNB8AxGRDhL1yY8jgPHAa2Y2Nyz7JVAK4O5TgO8A55pZA7AZGOvBIWsNZvYT4B9AMXCHu7+e7Q8Qta1bt3LTTTdx9NFHK6GISOQiTSru/jzJx0bi6/wJ+FOKadOB6R0QWt6ora3l/fffZ8qUKVGHIiIS/ZiKtM/111/PPvvsw/HHHx91KCIiSir5bNasWbz44ov89Kc/1aXtRSQnaEuUx2644QZ69+7NmWeeGXUoIiKAkkreWrFiBf/3f//HmWeeyY477hh1OCIigJJK3poyZQoNDQ389Kc/jToUEZHPKankoc8++4wpU6bwrW99i4EDB0YdjojI55RU8tB9993HypUrOe+8JncKEBGJlJJKnnF3rr/+evbff3+OOuqoqMMREdlO1GfUSys9//zzvPLKK0yZMoXgKjciIrlDeyp55vrrr2eXXXZh/PjxUYciItKEkkoeWbJkCQ8++CA/+MEPKCkpiTocEZEmlFTyyE033YSZ8eMf/zjqUEREklJSyRMbN27k1ltv5ZRTTtEdHUUkZymp5Im7776bNWvW6DBiEclpSip5wN254YYbGDJkCCNGjIg6HBGRlHRIcR546qmneOONN5g6daoOIxaRnKY9lTxw/fXXs/vuuzN27NioQxERSSvSpGJmA8yszswWmNnrZtZkwMDMKs1sXvj4l5kNipu2yMxeM7O5ZjY7u9Fnx9tvv82jjz7KOeecQ/fu3aMOR0Qkrai7vxqAi9z9ZTPrDcwxsyfd/Y24Ou8BR7r7GjM7DqgGhsdNH+XuH2Ux5qz605/+RNeuXTn33HOjDkVEpFlR36N+BbAifL3BzBYA/YE34ur8K67JDGCvrAYZofXr13PnnXdy+umns8cee0QdjohIs3JmTMXMyoFDgJfSVPs+8FjceweeMLM5ZlbVcdFF484772TDhg06jFhE8kbU3V8AmFkv4AHgfHdfn6LOKIKkckRc8Qh3X25muwNPmtmb7v5skrZVQBWQNycONjY2cuONN3L44YdTUVERdTgiIi3Soj0VM/tvMzu4IwIws64ECaXG3WtT1DkYuA04yd0/jpW7+/LweSXwIDAsWXt3r3b3Cnev6Nu3b6Y/QkbV1NRQXl5O165dWbhwIUOGDIk6JBGRFmtp99dU4OT4AjP7npk93Z6FW3DSxe3AAne/JkWdUqAWGO/ub8WV7xAO7mNmOwDfBOa3J56o1dTUUFVVxeLFi3F3AO644w5qamoijkxEpGXaM6ZSDhzZzuWPAMYD3wgPC55rZseb2Tlmdk5Y5zJgN+DPCYcO9wOeN7NXgZnAo+7+eDvjidTEiRPZtGnTdmWbNm1i4sSJEUUkItI6UR/99TyQ9hRxdz8bODtJ+bvAoKYt8teSJUtaVS4ikmty5ugvSX0QQb4cXCAioqSSQyZPnkyPHj22KyspKWHy5MkRRSQi0jqtSSreYVEIAJWVlZx00kkAmBllZWVUV1dTWVkZcWQiIi3TmjGVSWY2KbHQzBpT1Hd3z4nzYPLJe++9x9ChQ5k5c2bUoYiItFpr9lSslQ91rbXS0qVLmTlzJmPGjIk6FBGRNmnRnoS7K0FkwUMPPQSgpCIieUvJIofU1tZywAEH8KUvfSnqUERE2kRJJUesWrWKZ599VnspIpLXlFRyxLRp09i2bZuSiojkNSWVHFFbW8vee+/NoEEFdZEAEelklFRywLp163jqqacYM2YMwTU2RUTyk5JKDpg+fTpbtmxR15eI5D0llRxQW1vLHnvswVe/+tWoQxERaRcllYht3ryZ6dOnc8opp1BUpK9DRPKbtmIRe+KJJ9i0aZO6vkSkICipRKy2tpZddtmFI49s7/3ORESip6QSoa1btzJt2jROPPFEunbtGnU4IiLtFmlSMbMBZlZnZgvM7HUzOy9JHTOzG8zsHTObZ2ZD4qZ9z8zeDh/fy2707VdfX8/atWvV9SUiBSPqS9M3ABe5+8tm1huYY2ZPuvsbcXWOA/YNH8OBm4HhZrYr8GugguBeL3PMbJq7r8nuR2i72tpadthhB4455pioQxERyYhI91TcfYW7vxy+3gAsAPonVDsJ+IsHZgA7m9mewLHAk+6+OkwkTwKjsxh+uzQ2NvLggw9y/PHH07Nnz6jDERHJiJwZUzGzcuAQ4KWESf2BpXHvl4VlqcrzwowZM/jwww/V9SUiBSUnkoqZ9QIeAM539/WJk5M08TTlyeZfZWazzWz2qlWr2hdshtTW1tKtWzeOP/74qEMREcmYyJOKmXUlSCg17l6bpMoyYEDc+72A5WnKm3D3anevcPeKvn37ZibwdnB3amtrOeaYY9hxxx2jDkdEJGOiPvrLgNuBBe5+TYpq04D/Do8C+yqwzt1XAP8Avmlmu5jZLsA3w7KcN3fuXBYtWqSuLxEpOFEf/TUCGA+8ZmZzw7JfAqUA7j4FmA4cD7wDbALODKetNrMrgFlhu8vdfXUWY2+z2tpaioqKOPHEE6MORUQko8w96TBEwaqoqPDZs2dHGsMBBxxAv379ePrppyONQ0SkJcxsjrtXtKRu5GMqnc2bb77JG2+8oa4vESlISipZ9uCDDwJw8sknRxyJiEjmKalkWW1tLcOHD2evvfaKOhQRkYxTUsmiJUuWMHv2bHV9iUjBUlLJoljXl5KKiBQqJZUsqq2t5eCDD2bgwIFRhyIi0iGUVLLkww8/5LnnntNeiogUNCWVLJk2bRrurqQiIgVNSSVLamtrGThwIAceeGDUoYiIdBgllSxYu3Yt//znPxkzZgzB5c5ERAqTkkoWPPLII2zdulVdXyJS8JRUsqC2tpb+/fszdOjQqEMREelQSiodbOPGjTz++OOccsopFBVpdYtIYdNWrgPV1NRQXl7O5s2b+etf/0pNTU3UIYmIdKio76dSsGpqaqiqqmLTpk0ArFy5kqqqKgAqKyujDE1EpMNoT6WDTJw48fOEErNp0yYmTpwYUUQiIh1PSaWDLFmypFXlIiKFIOp71N9hZivNbH6K6b8ws7nhY76ZNZrZruG0RWb2Wjgt2ls5JlFaWtqqchGRQhD1nspUYHSqie5+tbsPdvfBwCXAMwn3oR8VTm/RbS6zafLkyU2O9iopKWHy5MkRRSQi0vEiTSru/iywutmKgXHAvR0YTkadeuqpmBm9e/fGzCgrK6O6ulqD9CJS0PLi6C8zKyHYo/lJXLEDT5iZA7e4e3UkwaXw4osv0tjYyD333MMJJ5wQdTgiIlmRF0kF+DbwQkLX1wh3X25muwNPmtmb4Z5PE2ZWBVRB9sY06uvrKSoq4mtf+1pWlicikguiHlNpqbEkdH25+/LweSXwIDAsVWN3r3b3Cnev6Nu3b4cGGlNXV8eQIUPYaaedsrI8EZFckPNJxcx2Ao4E/h5XtoOZ9Y69Br4JJD2CLAqbNm3ipZdeYtSoUVGHIiKSVZF2f5nZvcBIoI+ZLQN+DXQFcPcpYbVTgCfcfWNc037Ag+Fl5LsA97j749mKuzkvvvgiW7ZsUVIRkU4n0qTi7uNaUGcqwaHH8WXvAoM6Jqr2q6uro7i4mCOOOCLqUEREsirnu7/yUX19PRUVFfTu3TvqUEREskpJJcM2btzIzJkz1fUlIp2SkkqGvfDCC2zdupWRI0dGHYqISNYpqWRYfX09Xbp0YcSIEVGHIiKSdUoqGVZXV8ewYcPo1atX1KGIiGSdkkoGbdiwgVmzZqnrS0Q6LSWVDHrhhRdobGzUIL2IdFpKKhlUV1dH165dOfzww6MORUQkEkoqGVRXV8fw4cMpKSmJOhQRkUgoqWTI+vXrmTNnjrq+RKRTU1LJkOeee45t27ZpkF5EOjUllQypr6+nW7duHHbYYVGHIiISGSWVDKmrq+Owww6jZ8+eUYciIhIZJZUMWLt2La+88oq6vkSk01NSyYDYeIoG6UWks1NSyYC6ujp69OjB8OHDow5FRCRSSioZEBtP6dGjR9ShiIhEKtKkYmZ3mNlKM0t6f3kzG2lm68xsbvi4LG7aaDP7t5m9Y2YTshf19lavXs2rr76qri8REaLfU5kKjG6mznPuPjh8XA5gZsXATcBxwP7AODPbv0MjTeHZZ5/F3ZVURESIOKm4+7PA6jY0HQa84+7vuvsW4D7gpIwG10J1dXX07NmToUOHRrF4EZGcEvWeSkscZmavmtljZnZAWNYfWBpXZ1lYlnX19fWMGDGC7t27R7F4EZGckutJ5WWgzN0HATcCD4XllqSup5qJmVWZ2Wwzm71q1aqMBffRRx8xb948dX2JiIRyOqm4+3p3/yR8PR3oamZ9CPZMBsRV3QtYnmY+1e5e4e4Vffv2zVh8zzzzDIBOehQRCeV0UjGzPczMwtfDCOL9GJgF7Gtme5tZN2AsMC3b8dXX11NSUqLxFBGRUJcoF25m9wIjgT5mtgz4NdAVwN2nAN8BzjWzBmAzMNbdHWgws58A/wCKgTvc/fVsx19XV8cRRxxB165ds71oEZGcFGlScfdxzUz/E/CnFNOmA9M7Iq6WWLlyJa+//jrf/e53owpBRCTn5HT3Vy6LjadokF5E5AtKKm1UV1dHr169GDJkSNShiIjkDCWVNqqvr+drX/uaxlNEROIoqbTBBx98wIIFC9T1JSKSQEmlDerr6wGdnyIikkhJpQ3q6+vZcccdOeSQQ6IORUQkpyiptEFdXR1f//rX6dIl0iOyRURyjpJKKy1fvpy33npLXV8iIkkoqbRSbDxFg/QiIk0pqbRSXV0dO++8M4MGDYo6FBGRnKOk0kqx8ZTi4uKoQxERyTlKKq2wdOlSFi5cqK4vEZEUlFRaQeeniIikp6TSCvX19ey6664cfPDBUYciIpKTlFRaoa6ujiOPPJKiIq02EZFktHVsocWLF/Pee++p60tEJA0llRaoqanh0EMPBeCqq66ipqYm4ohERHJTpEnFzO4ws5VmNj/F9Eozmxc+/mVmg+KmLTKz18xsrpnN7qgYa2pqqKqq4uOPPwZgxYoVVFVVKbGIiCQR9Z7KVGB0munvAUe6+8HAFUB1wvRR7j7Y3Ss6KD4mTpzIpk2btivbtGkTEydO7KhFiojkrajvUf+smZWnmf6vuLczgL06OqZES5YsaVW5iEhnFvWeSmt8H3jAjrP9AAAO50lEQVQs7r0DT5jZHDOr6qiFlpaWtqpcRKQzy4ukYmajCJLK/8QVj3D3IcBxwI/N7Otp2leZ2Wwzm71q1apWLXvy5MmUlJRsV1ZSUsLkyZNbNR8Rkc4g55OKmR0M3Aac5O4fx8rdfXn4vBJ4EBiWah7uXu3uFe5e0bdv31Ytv7KykurqasrKyjAzysrKqK6uprKysk2fR0SkkOX0XabMrBSoBca7+1tx5TsARe6+IXz9TeDyjoqjsrJSSUREpAUiTSpmdi8wEuhjZsuAXwNdAdx9CnAZsBvwZzMDaAiP9OoHPBiWdQHucffHs/4BRERkO1Ef/TWumelnA2cnKX8X0A1NRERyTM6PqYiISP5QUhERkYxRUhERkYwxd486hqwys1XA4jY27wN8lMFwCoHWSVNaJ01pnTSVT+ukzN1bdD5Gp0sq7WFmszvyOmP5SOukKa2TprROmirUdaLuLxERyRglFRERyRglldZJvPS+aJ0ko3XSlNZJUwW5TjSmIiIiGaM9FRERyRgllRYws9Fm9m8ze8fMJkQdT67I1i2dc1myW2Kb2a5m9qSZvR0+7xJljNmWYp1MMrP3w9/KXDM7PsoYs83MBphZnZktMLPXzey8sLzgfitKKs0ws2LgJoL7tuwPjDOz/aONKqd0+C2dc9xUmt4SewLwT3ffF/hn+L4zmUry24RfG/5WBrv79CzHFLUG4CJ3/wrwVYJ7QO1PAf5WlFSaNwx4x93fdfctwH3ASRHHJDnC3Z8FVicUnwTcFb6+Czg5q0FFLMU66dTcfYW7vxy+3gAsAPpTgL8VJZXm9QeWxr1fFpZJlm7pnIf6ufsKCDYmwO4Rx5MrfmJm88Lusbzv5mkrMysHDgFeogB/K0oqzbMkZTpkLtDiWzpLp3czsA8wGFgB/DHacKJhZr2AB4Dz3X191PF0BCWV5i0DBsS93wtYHlEsOaU1t3TuZD40sz0BwueVEccTOXf/0N0b3X0bcCud8LdiZl0JEkqNu9eGxQX3W1FSad4sYF8z29vMugFjgWkRxxQ5M9vBzHrHXhPc0nl++ladxjTge+Hr7wF/jzCWnBDbcIZOoZP9Viy4Te3twAJ3vyZuUsH9VnTyYwuEhz9eBxQDd7j75IhDipyZ/SfB3gl8cUvnTrde4m+JDXxIcEvsh4C/AqXAEuC/3L3TDFynWCcjCbq+HFgE/DA2ltAZmNkRwHPAa8C2sPiXBOMqBfVbUVIREZGMUfeXiIhkjJKKiIhkjJKKiIhkjJKKiIhkjJKKiIhkjJKK5BwzqzczN7NJUccSJTMrNrMLzewVM9sYrhM3sxZdH0rr8QtmVh63/sqjjqeQKankifDS4bE/io1m9h9p6sb/AY3MYpiSWdcRXM5kMMG5QB+Gj0/bO2MzOz/8TQ1u77yiFn6OSUoWuaFL1AFIm5QQnFD2w6gDkY4RXq0g9v1eDPzBW39S2RLg38BHSaadD5QRnIg4t41h5opfh8/1BJ8nma0E6yL2WjqIkkr+OsvM/ujub0UdiHSILwNdw9c3tyGh4O7/ndmQ8pe7v0+wTqWDqfsr/ywF5hH8Q3BlxLFIxymJvXD3T6IMRKQ1lFTyzzbgkvD1qWbWqqu9tnTAMrxVsJvZGenam1mZmd1qZkvM7FMzW2hmvw0vMhlrc6CZ3W1mS8M6b5vZr8KrtjYXbzczmxDeh2Ojma0Jb7t6XAva7mNmN4a3cP3EzDaFr68zs9IUbc4IP9ui8P0oM3vIzFaYWaOZTW1uuQnzKzazs8zsaTP7yMw+s+C2uvcnG++KLZ+gKydW5nGP+sQ2aZbdZKA+NjZH0PUFcGfC/JPuEZnZSDO7N+57XmdmM83s4vjvOqHN1HCeUy1wtpk9b2YfJ/62zGyImV1mZs+a2eJwGWvNbIaZ/Y8Fl4xPOv+4orqEz7Iorm6zv3sz2ymM4WUzW29mm8Pf6s0WXOsu1Xr+fPzSzHqHv/83w/Yfm9kjZjY8VfuC4+565MEDmER4Mb7wfX34/ukkdcvDaQ6MTDOtPM3yFoV1zkjTfgywJny9juCWqbFpzxJ033wL2BiWrSVIirE696VYduyzXRnOxwn6wdfEtXVgUpr4fwBsiav7KbAp7v064Jgk7c7gi4se/iwu3rXh/Ka24jvbCaiLW2ZD+Bni18HVCW1OBz4guHNirM4HcY/aViw/th4nxZX9PJxPY9x6iJ//Bwnz6EJwqfr49b4h4bt+EyhLsvyp4fS7gPvD143hZ2uM/20lzL8xyXf9OrB7wvyvD2OO1Vmd8FlmtfR3DxxA0AsQq7MZWJ/w+zk1xXqO1RkHvB3XfmPctC3AsVFvR7LxiDwAPVr4RTVNKsPjfrCjE+rG/wGNTDOtPM3yFtF8UlkDPAXsH07rCfw0boNzBcHG+L7YRgfoBfw2bh5HJ1l2PV9syD8lGLDuEU4bELeBcuDEJO1PjvtD/l+C/8otfOxHcFXY2Aa1NKHtGXEbhQbgTmBAOK0Y2KcV39nfwnl9Fq6XkrB8D4LLoMc+wzlJ2o6MTW/Hbya2Hie19PtNUu86vkhs5wK7huVdwxhfDqfPAYoS2k7liyS0FbgI2DHud7BnXN0ngTMJrtbbJe73dApB0nJSJFRS/NZb+rsHegPvhtOWAcfHPgswCHiRLxLLoDTLX02Q/EYR9AIZMDQu/kWJ66gQH5EHoEcLv6iEpBKW1YZlrxBecTosj/8DGpkwn5R/XAn1km50EtrPB7onafuXuDpPxMcWVye2B3Jbkmn1ce3PSjK9CHgmnP56wrRu4YYhadu4en8P61yXUH5G3LIfaMf3NSxuPlUp6sSSzirCpBk3bWSsfTtiiK3HSS39fhPqHEiwV7UROChFnd588R/+yQnTpsatg5+243P0J9igbyPhn4BwenuTyv/wxT8hB6b4jO+FdR5Js/yVJOxNhdMPiqszoq3rIV8eGlPJb78k6CoYTLDrnW3XuvtnScr/Eff6Kg//slLUOTjN/JcS7Clsx4O7B/42fLu/mR0UN/k4go3Qh8naxvlL+Hxsmjr/m2Zac8aGz8uA21LUuTR87gMc045ldZTvE/y3/ai7v5asgrtvILh/DKRel2uAW9oahAdHbr0axnJ4W+eTxunh89/cvcnNw8LP+Pvw7XFmtlOK+VR7cBfUxPavESQlSP97Lwg6pDiPufubZnYncDZwhZnd7+7ZPAZ/ZoryD+Nez2qmzi5p5l+fIiFBsKfTQPAbriC4+RHAEXHzXWFmqebdLXwuSzF9M0HXTltVhM91YRJswt0XmNn7BEmwAni4HcvrCLF1eZyZfZCmXmwQPdW6nOXuW9ItyMyKCBLxWIJ/kvoCPZJU3SvdfFrLgru5xjb0T6Wp+mT4XAQMIRgrS/RSmvbLgb2BXVsbY75RUsl/k4BK4D+Bc4Abs7jsDSnKG2Ivwv/y0tVJdwTY+6kmuPtnZvYx0A/YPW5S7EoD3cJpzemZovzjVMmghWIxpfwMoWUESWX3ZupFIbYue/FF4kinJEV52vuum1kJ8AjBWETMFoIxitg/SbsS/FaSHmnWDrsSjJVB+u9qWdzrVN9Vqt86tOz3XhDU/ZXnwq6BWCL5VbJDL/NYqr2UdGIbiMfd3VrySDGfxjbGnKiln6Etn7WjxdblhBauy5Ep5tPcupxIkFA2AxcQ7PH0cPfd3H0Pd9+DL/YCUu56ZkC678BTvJYESiqF4X8J+q13JzjCJp2GuNfJuhdiUvUbZ1PKrg4z6w7sFr6N/0841k1zENGKxTSgmXqxz7iqA2Npq2yty9j40+Xufp27L0nS7blHBy07dngzpP+u4qfl4neVM5RUCoC7rwWuCt9eRPqulDVxr5P+EZnZl4CdMxNduxxpqQdFvsYX3bez48pfCJ/7m9kRRCcW06hwvKAJM/syQdcXpB576iixrr10//nH1uW3OngPOPY7fCXZxPBkxYFp2scSUKv3YsKxnnnh26PSVD06fN5G+8baCp6SSuG4gaDftzfwq1SV3H0jsDB8e2qKahMzG1qblQLfSywMN9K/DN8uSDgy6WFgRfj6+rC/PiUz66iB0/vC5/4EB1Ikc3n4/BHpB4k7wvrwOd0/D7GTHncGrk43MzPr2o7Esy58HpRi+lUpymNa8lnSiX1X3zGzAxMnhp/r4vDtdHdfl1hHvqCkUiDc/VOCQXuAbzdT/d7w+Swz+5GZ9QQwswFmdhvBIZabOiTQ1lkH3GxmPzCzHhDESBB/bFB3uwQYrocfEWwMhwAvmNmx4VE+hPPY28x+aGYzw7oZ5+4zgQfCtzea2U9iCc7M9jCzW4H/CqdfGsadTbFDZ79jZkmPwHP3uQQnPwKcE15aZnBs79GCS9AMMrNLCf5Raetl9B8Pn39lZmPMrEs4/73N7B7gNLbfw071WSqb+ycihZsJDvntCjxmZsfF9i7Dw9X/QXDk1hbS/MMmoahPlNGjZQ+SnPyYpE4xsIDtL28xMkm9XgRn/sbqxF8WYwtBH/cimj/5sTxFHCNjddLEekaqz8P2l2l5Li6u+EuXOHBFmvlXsv1lMrYS7BF8mjCPiS2Nqw3f2U5sfyLn1vAzpLxMS2vWYQuWH1v2pCTTvh4XRwPBIa+LEj93+Ju6NmGdbQ7X5daE8hEJbaeG5VObibOM7S+3spXgagqx95c081m+G1d3C8Ee+yLg+Zb+bglO9FzG9p9xXdz7T4HvpIi/JSdfpoy/0B7aUykg7t7IF91C6ep9QnAOwjUE/6E1EPwhPwAc5u73pWmeTVsI+rl/SXAvjO4Ef+j/BL7l7pemaujuNQT98L8lGN/4hKB75FOC+4f8iaCf/HcdFbwH3SRHEZxEWE9wyGkvgg3oA8Aod/9FRy2/mdieJbgu21ME67Qfwca9LKFeo7tfQLDXV03wPTQSJMw1BOMuk4DB7v4CbeDuiwnO07mdILlB8D09QnC9rLQnobr73cB44HmCPew9w8/R4nNaPDjp8QCCzzKX4G+iO8Ee2BTgAHf/W4s/VCdmYRYVERFpN+2piIhIxiipiIhIxiipiIhIxiipiIhIxiipiIhIxiipiIhIxiipiIhIxiipiIhIxiipiIhIxiipiIhIxiipiIhIxvw/YMCWsR3dJOcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x2c92c3b3128>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.plot(F,'-ok')\n",
    "plt.xlabel('Number of iteration',fontsize=25)\n",
    "plt.ylabel('F',fontsize=20)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "图4.9显示了$F$随着时间的变化。在大约12次迭代之后，$F$变得几乎恒定。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
